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ABSTRACT 

Large scale structure deflects cosmic microwave background (CMB) photons. Since large angular 
scales in the large scale structure contribute significantly to the gravitational lensing effect, a realistic 
simulation of CMB lensing requires a sufficiently large sky area. We describe simulations that include 
these effects, and present both effective and multiple plane ray-tracing versions of the algorithm, 
which employs spherical harmonic space and does not use the flat sky approximation. We simulate 
lensed CMB maps with an angular resolution of ~ 0/9. The angular power spectrum of the simulated 
sky agrees well with analytical predictions. Maps generated in this manner are a useful tool for the 
analysis and interpretation of upcoming CMB experiments such as PLANCK and ACT. 
Subject headings: cosmic microwave background — gravitational lensing — large-scale structure of 
universe — methods: N-body simulations — methods: numerical 
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1. INTRODUCTION 

While the current generation of CMB experiments have 
had a significant impact on cosmology by helping to es- 
tablish a standard paradigm for cosmology (Spergel et al. 
2003, 2007), the upcoming generation of CMB experi- 
ments still has the potential to provide novel new in- 
sights into cosmology. PLANCK 1 and ground based ex- 
periments, such as the Atacama Cosmology Telescope 
(ACT) 2 , will be mapping the CMB sky with significantly 
higher angular resolution than ever before. Secondary 
anisotropies on small angular scales encode important in- 
formation about the late time interaction of CMB pho- 
tons with structure in the Universe. One of the most 
basic of these interactions is the gravitational effect of 
the large scale structure potentials deflecting the paths 
of the photons, an effect justifiably referred to as the 
Gravitational Lensing of the CMB. 

The effect of gravitational lensing can be thought of as 
a remapping of the unlensed CMB field by a line-of-sight 
averaged deflection field (for a recent review, see Lewis 
& Challinor 2006). Therefore, lensing does not change 
the one-point properties of the CMB. However, it does 
modify the two and higher-point statistics, and gener- 
ates non-Gaussianity (Seljak 1996; Zaldarriaga & Seljak 
1999; Zaldarriaga 2000). Although the typical deflection 
suffered by a CMB photon during its cosmic journey is 
about three arcminutes, the deflections themselves are 
coherent over several degrees, which is comparable to the 
typical size of the acoustic features on the CMB. Thus 
lensing causes coherent distortions of the hot and cold 
spots on the CMB, and thereby broadens their size dis- 
tribution. This leads to redistribution of power among 
the acoustic scales in the CMB, and shows up in the 
two-point statistics as a smoothing of the acoustic peaks. 
At smaller scales, where the primordial CMB is well ap- 
proximated by a local gradient, deflectors of small an- 
gular size produce small-scale distortions in the CMB, 

Electronic address: sudeep@astro.princeton.edu,bode@astro.princeton 

1 http:/ /www.rssd.esa.int/index.php?project=Planck 

2 http:/ /www. physics. princeton.edu/act 



thereby transferring power from large scales in the CMB 
to the higher multipoles. Also, although the primordial 
CMB can be safely assumed to be a Gaussian random 
field (Komatsu et al. 2003), and the large scale lensing 
potential can also be well approximated by a Gaussian 
random field, the lensed CMB — being a reprocessing 
of one Gaussian random field by another — is itself not 
Gaussian. The effect of lensing on the power spectrum 
of the CMB is important enough that it should be taken 
into account while deriving parameter constraints with 
future higher resolution experiments. But what is even 
more interesting is that the non-Gaussianity in the lensed 
CMB field should enable us to extract information about 
the projected large scale structure potential, and thereby 
constrain the late time evolution of the Universe and 
Dark Energy properties. Therein lies the main moti- 
vation of studying this effect in utmost detail. Progress 
in this area has been slow. Measurements of the CMB 
precise enough to enable a detection of weak lensing were 
not available in the pre-WMAP era. Also, picking out 
non-Gaussian signatures in the measured CMB sky by 
itself is extremely difficult, due to confusion from sys- 
tematics, foregrounds, and limited angular resolution. 

Rather than looking at signatures of lensing only in the 
CMB, one can also measure to what extent the deflection 
field estimated from the CMB correlates with tracers of 
the large scale structure which contributed to the lensing. 
It is easily realized that this approach is powerful (Peiris 
& Spergel 2000) because many of the systematics disap- 
pear upon cross-correlating data sets. This approach was 
taken in recent years by Hirata et al. (2004) and Smith 
et al. (2007), using WMAP 1-year and 3-year data re- 
spectively. The former work looked at the cross correla- 
tion with SDSS luminous red galaxies (LRG), while the 
latter used the NRAO-VLA Sky Survey (NVSS) radio 
sources as their large scale structure tracers. As the lens- 
ing efficiency for the CMB is highest between redshifts of 
one and four, higher redshift tracers should show greater 
edsross correlation signal, which makes the NVSS radio 
sources better tracers for such study; Smith et al. (2007) 
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report a 3.4<r detection. An independent analysis by Hi- 
rata et al. (2008) looking for this effect in the WMAP 
3-year data in cross correlation with SDSS LRG+QSO 
and NVSS sources find this signal at the 2.5a level. With 
these pioneering efforts and with higher resolution CMB 
data from experiments such as ACT, PLANCK and the 
South Pole Telescope (SPT) 3 on the horizon, we are en- 
tering an era where robust detection and characterization 
of this effect will become a reality. Also, with upcom- 
ing and proposed large scale structure projects (LSST 4 , 
SNAP 5 , ADEPT 6 , DESTINY 7 , etc.) there will in future 
be many more datasets to cross-correlate with the CMB. 

One of the immediate results of such cross-correlation 
studies will be a measurement of the bias of the tracer 
population. Because such cross correlations tie together 
early universe physics from the CMB and late time evolu- 
tion from large scale structure, they will also be sensitive 
to Dark Energy parameters (Hu et al. 2006) and neutrino 
properties (Smith et al. 2006a; Lesgourgues et al. 2006), 
and can potentially break several parameter degenera- 
cies in the primordial CMB (M. Santolini, S. Das and 
D. N. Spcrgcl, in preparation). Combination of galaxy 
or cluster lensing of the CMB with shear measurements 
from weak lensing of galaxies can also provide impor- 
tant constraints on the geometry of the Universe (S. Das 
and D. N. Spergel, in prep; Hu et al. 2007). Again, 
with high enough precision of CMB data, it is possi- 
ble to estimate, using quadratic (Okamoto & Hu 2003) 
or maximum likelihood (Hirata & Seljak 2003) estima- 
tors, the deflection field that caused the lensing. Such 
estimates can be turned into strong constraints of the 
power spectrum of the projected lensing potential (Hu 
& Okamoto 2002), which is also sensitive to the details 
of growth of structure. The estimated potential from 
the lensed CMB alone, or the potential estimated from 
weak lensing surveys (Marian & Bernstein 2007) , can be 
also used to significantly de-lens the CMB. This is par- 
ticularly important in the detection of primordial tensor 
modes via measurements of CMB polarization. This is 
because (even though detection of the so-called B modes 
in CMB polarization is hailed as the definitive indicator 
of the presence of gravitational waves from the inflation- 
ary era) these mode can be potentially contaminated by 
the conversion of E-modes into B-modes via gravitational 
lensing. De-lensing provides a way of cleaning these con- 
taminating B-modes produced by lensing and thereby 
probing the true gravitational wave signature. 

In this paper, we describe a method for simulating the 
gravitational lensing of the CMB temperature field on a 
large area of the sky using a high resolution Tree-Particle- 
Mesh (TPM; Bode et al. 2000; Bode & Ostriker 2003) 
simulation of large scale structure to produce the lens- 
ing potential. The reason for considering a large area of 
the sky is twofold. First, the deflection field has most 
of its power on large scales (the power spectrum of the 
deflection field peaks at I ~ 50 in the best-fit cosmologi- 
cal model), and much of the power redistribution in the 
acoustic peaks of the CMB occurs via coupling of modes 

3 http://spt. uchicago . edu 

4 http:/ /www. lsst.org/lsst_homc.shtml 

5 http://snap.lbl.gov/ 

6 http:/ /universe. nasa.gov/program/probcs/adcpt. html 

7 http://destiny.asu.edu/ 



in the CMB with these large coherent modes in the de- 
flection field. A large sky allows for several such modes 
to be realized. It is estimated that a small (flat) sky sim- 
ulation that misses these modes would typically under- 
estimate the lensing effect by about 10% in the acoustic 
regime, and more in the damping tail (Hu 2000). Sec- 
ond, one of the major goals of simulations such as this is 
to produce mock observations for upcoming CMB exper- 
iments. PLANCK is an all-sky experiment, and many of 
the future CMB experiments (including ACT and SPT) 
will observe relatively large patches of the sky. There- 
fore, simulating CMB fields on a large area of the sky 
is a necessity. This method fully takes into account the 
curvature of the sky. Although presented here for a polar 
cap like area, it can be trivially extended to the full sky. 

The value of a simulation as described here is multi- 
faceted, particularly in the development of algorithms for 
detection and characterization of the CMB lensing effect 
for a specific experiment. Since each experiment has a 
unique scanning mode, beam pattern, area coverage, and 
foregrounds, operations and optimizations performed on 
the data to extract the lensing information will have to 
be tailored to the specific experiment. A large-sky lensed 
CMB map acting as an input for a telescope simula- 
tor provides the flexibility of exploring various observing 
strategies, and also allows for superposition of known 
foregrounds. Another important aspect of this simula- 
tion is that the halos identified in the large scale struc- 
ture simulation can be populated with different tracers 
of interest. Also, other signals, such as the Thermal and 
Kinetic Sunyaev Zel'dovich effects and weak lensing of 
galaxies by large scale structure, can be simulated using 
the same large scale structure. This opens up the possi- 
bility of studying the cross-correlation of the CMB lens- 
ing signal with various indicators of mass, and thereby 
predicting the level of scientific impact that a specific 
combination of experiments can have. 

As noted in Lewis (2005), the exact simulation of the 
lensed CMB sky, which requires the computation of spin 
spherical harmonics on an irregular grid defined by the 
original positions of the photons on the CMB surface, 
is computationally expensive and requires robust par- 
allclization. Lewis (2005) suggested an alternative in 
which one would rcsamplc an unlenscd CMB sky, gen- 
erated with finer pixelation, at these unlensed positions. 
This method was implemented in the publicly available 
LcnsPix 8 code that was based on Lewis (2005). How- 
ever, producing a high resolution lensed map requires a 
much higher resolution unlensed map, the generation of 
which becomes computationally more expensive as reso- 
lution increases. Here we put forward another alterna- 
tive, in which we do the resampling with a combination 
of fast spherical harmonic transform on a regular grid 
followed by a high order polynomial interpolation. This 
interpolation scheme has been adapted from Hirata et al. 
(2004), and is called the Non-Isolatitude Spherical Har- 
monic Transform (NISHT). This method is accurate as 
well as fast, and does not require parallclization or pro- 
duction of maps at a higher resolution. Another added 
advantage of this method is that the same algorithm can 
be used to generate the gradient of a scalar field on an 
irregular grid. Since the deflection field is a gradient 

8 http://cosmologist.info/lcnspix 
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of the lensing potential, this opens up the possibility of 
performing a multiple plane ray tracing simulation. This 
is because the rays, as they propagate from one plane to 
another, end up on irregular grids, so the deflection fields 
on the subsequent planes have to be evaluated on irregu- 
lar grids. At the time of the development of this project, 
LensPix did not include an interpolation scheme, and 
used the methods as described originally in that paper. 
Concurrently with the completing of the current work, 
an interpolation scheme (Akima 1996) different from the 
one described here has been added to that code. An- 
other notable difference of our results with LensPix, is 
that while the latter uses a Gaussian Random realization 
of the deflection field, we have used a large scale struc- 
ture simulation to produce the same, thereby including 
all higher order correlations due to non-linearities. 

The paper is laid out as follows. In §2 we explain the 
lensing algorithm, describing the governing equations in 
§2.1 and their discretization in §2.2. Then we discuss the 
effective lensing approach (§2.3) as well as the multiple 
plane ray tracing approach (§2.4). At the heart of the 
lensing algorithm lies the non-isolatitude spherical har- 
monic transform algorithm adapted from Hirata et al. 
(2004), which is reproduced in some detail for complete- 
ness in §2.5. As discussed earlier, we have employed a 
light cone A-body simulation and adopted a special po- 
lar cap like geometry for generating the lensing planes 
(§3). For comparison of the simulated fields with theo- 
retical prediction, we compute the angular power spectra 
on the polar cap window; in §4 we describe some of the 
subtleties involved in computing the power spectra. We 
present our results in §5 and describe the tests that we 
have performed in §6. Conclusions are presented in §7. 

2. THE LENSING ALGORITHM 

2.1. Basic Equations 

We would like to note here that while the calculations 
for the simulation described here has been done for a 
flat universe, our approach is generalizable to non-flat 
geometries. 

The deflection angle of a light ray propagating through 
the space is 

da = -2V±_^dn, (1) 

where is da is the deflection angle, IP is the Newtonian 
potential, Vj_ denotes the spatial gradient on a plane 
perpendicular to light propagation direction and n is the 
radial comoving distance. The transverse shift of the 
light ray position at n due to a deflection at rj is given 

by 

d X (n)=d A (n-n')da(n'), (2) 

where dA{rj) is the comoving angular diameter distance. 

The final angular position 6{rj) = x(r]) / dA{rj) is there- 
fore given by 

= 9(0) + £.(>,), (3) 
where a. is the total effective deflection. 

2.2. Discretization 

We will now discretize the above equations by dividing 
the radial interval between the observer and the source 



/ dr]'d A (ri-ri')V±y 
Jo 



into N concentric shells each of comoving thickness A77. 
We project the matter in the i-th shell onto a spherical 
sheet at comoving distance rji which is halfway between 
the the edges of the shell (i increases as one moves away 
from the observer) . Since we shall be working in spherical 
coordinates it is advantageous to use angular differential 
operators instead of spatial ones. We rewrite equation 
(1) in terms of the angular gradient Vj as 



da 



\7 h ^dn. 



(4) 



At the j-th shell at rjj, the deflection angle due to the 
matter in the shell can be approximated by an integral 
of the above: 



cr 



-f 



■ nj +A n /2 



= -V fi ^'(n), 



A77/2 



W & ^(fjh;fj)dfj 



(5) 
(6) 



where we have defined the 2-D potential on the sphere 

as 

= Tl — T / ^{fjn;fj)dfj. (7) 

d A{Vj) Jrn-Ari/2 

Here, the notation (nh; rj) signifies that the potential is 
evaluated at the conformal look-back time n 1 when the 
photon was at the position nh. The potential can be 
related to the mass overdensity in the shell via Poisson's 
equation, which reads 



V 2 * = 



AttG p - p 
c 2 (1 + zf 



(8) 



p being the mean matter density of the universe at red- 
shift z. By integrating the above equation along the line 
of sight, one can arrive at a two dimensional version of 
the Poisson equation (Vale & White 2003), 



where the surface mass density 



Vj+A-n/2 



(p - p)di). 



(9) 



(10) 



r )j -Ar ? /2 



Note that in going from the three dimensional to the 
two dimensional version, the term containing the radial 
derivatives of the Laplacian can be neglected (Jain et al. 
2000) . One can show that this term is small by expanding 
the potential W in Fourier modes k, with components fcii 
parallel to the line of sight and k± transverse to it. Then, 
the ratio of the components of the line of sight integral 
in the parallel and transverse directions will be ~ fc 2 /fcj_- 
Due to cancellation along the line of sight, only the modes 
with wavelengths comparable to the line of sight depth 
of each slice will survive the radial integral. These would 
be the modes with k\\ < On the other hand, the 

transverse component gets most of its contribution from 
scales smaller than ~ 100 Mpc i.e. k± ^> 27r/100 <~ 
0.1 Mpc -1 . Under the effective lensing approximation, 
the projection is along the entire line of sight from zero 
rcdshift to the last scattering surface, An ~ 10 4 Mpc, 
giving km < 10~ 3 Mpc -1 . Therefore, in this case the ratio 
of the radial and transverse components of the integral 
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Fig. 1. — Geometry illustrating the point remapping used in the 
text 

will be ~ k?,/kj_ "C 10~ 4 . For a multiple plane case, we 
would typically employ 10 lensing planes for which this 
ratio would be <C 1CP 2 . The approximation will break 
down if we employ thin shells. 
Defining the field K as 



4irG d,A(r]j 



C 2 (i+ Zj r 1 



A J E (n), 



(11) 



equation (9) takes the form 

V|^'(n) = 2K i (n). 

It is convenient to define an angular surface mass density 
A|(n) as the mass per steradian, 



(12) 



f 

A«?(n) = / 



(1 + z) 3 



(13) 



The surface mass density defined in equation (10) is re- 
lated to this through the relation 



This implies the following form of equation (11), 



AttG (1 + Zj) 
c 2 d A {r]j) 



(14) 



(15) 



Equation (15) is the key equation here. The quantity 
K can be readily calculated once the mass density is 
radially projected onto the spherical sheet. Expanding 
both sides of equation (12) in spherical harmonics, one 
has the following relation between the components: 



1(1 + 1) 



(16) 



It is interesting to note that the apparently divergent 
monopole (/ = 0) modes in the lensing potential can be 
safely set to zero in all calculations, because a monopole 
term in the lensing potential does not contribute to the 
deflection field. Being the transverse gradient of the po- 
tential, the deflection angle a(n) is a vector (spin 1) 
field defined on the sphere and can be synthesized from 
the spherical harmonic components of the potential in 
terms of vector spherical harmonics, as will be described 
in § 2.5. 

2.3. Connection with effective lensing quantities 

In weak lensing calculations, one often takes an effec- 
tive approach, in which one approximates the effect of 



deflectors along the entire line of sight by a projected 
potential or a convergence which is computed along a 
fiducial undeflected ray (often referred to as the Born 
approximation). One therefore defines an effective lens- 
ing potential out to comoving distance n s as 



(17) 



In terms of the projected potential, the effective de- 
flection (see Eq. 3) is given by the angular gradient, 
a = -Vft^^. An effective convergence is also defined 
in a similar manner: 

«(») = ^V|0 e// (n) 



2 

=/*, 



^.-yv ^;,,, (18) 

In terms of the fields and K 3 defined on the multiple 
planes, these quantities are immediately identified as the 
following sums, 



k e f f 



(A)*E 



dAjh - Vj) 



^(n), 



(19) 
(20) 



Once k is obtained one can go through the analog of 
equation (16) and take its transverse gradient to obtain 
the effective deflection a. Using this, one can find the 
source position corresponding to the observed position 
0(0): 

6 s = 6{0) + a. (21) 

In §5, we shall use this effective or single plane approxi- 
mation to lens the CMB. 

Equation (21) is to be interpreted in the following man- 
ner (Challinor & Chon 2002). The effective deflection 
angle is a tangent vector at the undeflected position of 
the ray. The original position of the ray on the source, or 
unlcnscd, plane is to be found by moving along a geodesic 
on the sphere in the direction of the tangent vector and 
covering a length a of an arc. The correct remapping 
equations can be easily derived from identities of spheri- 
cal triangles (Lewis 2005). For completeness, we give the 
derivation here. 

In Fig. 1, let the initial and final position of the ray in 
question be the points A= (9,4>) and B= (8',<j> + A</>), 
respectively. The North pole of the sphere is indicated 
as C, so that the dihedral angle at A is also the angle 
between a and —eg, so that 



a = —a cos See + & sin Jew,. 



(22) 



Now, applying the spherical cosine rule to the triangle 
ABC, we have 

cos 9' = cos 9 cos a + sin 9 sin a cos 6, (23) 

and applying the sine rule 



• a , . _ sin<J 
sin Aa> = sin a— — — . 

sin 0' 



(24) 



We use these equations to remap points on the CMB sky 
and on the intermediate spherical shells in the multiple 
plane described below. 
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Fig. 2. 
method. 



Geometry illustrating the multiple plane ray tracing 



2.4. Multiple plane ray tracing 

In the multiple plane case, we shoot ray outwards from 
the common center of the spherical shell (i.e. the ob- 
server) and follow their trajectories out to the CMB 
plane, thereby studying the time reversed version of the 
actual phenomenon. We assume all intermediate deflec- 
tions are small, as is really the case. Here we describe 
how we keep track of a ray propagating between multiple 
planes, as shown in Fig. 2. We assume a flat cosmology 
for this purpose. At some intermediate stage of the ray 
propagation, let a ray be incident on the i-th plane at the 
point A, where it gets deflected and reaches the i+l-th 
plane at the point D. The ray incident at A will not in 
general lie on the same plane as defined by the deflected 
ray AD and the center O of the sphere, which wc also 
consider as the plane of the figure. Assuming that we 
know the incidence angle /3i, we can obtain the addi- 
tional angle of deflection on due to the matter on plane 
i and compute the net deflection ot{ + (3^ Let us de- 
note by on, the effective angle, by which a ray has to be 
remapped from its observed position 6(0) to its current 
position 8i on plane i, so that 



(25) 



Obviously, di — and Q\ — 9(0). Therefore, the ef- 
fective angle (a,+i — a,) by which the ray has to be 
remapped from point B to point D on the shell i+l can 
be readily calculated from two descriptions of the arc 
BD, 



r?i+i(ai+i - a,) = (r]i+i -%)(<*,+ A). 



(26) 



In order to repeat this process for the (i+2)-th shell, one 
needs to know the value of the new incidence angle 
We now equate two ways of finding the length of the arc 
AC, 

77,(0;,+! - dj) = (n i+ i - r/t)Pi+i- (27) 
Substituting (d i+ i — dj) from equation (26), 

Vi 



Vi+i 



(28) 



Since we shoot the rays radially on the first plane, 
f3 1 = 0; therefore equations (26) and (28) can be used 
to propagate the ray back to the CMB surface, which 



we take to be the (N + l)-th plane, i.e. S = 0n+i- 
Although we only discuss results obtained with the ef- 
fective or single plane approximation here, the multiple 
plane version is straightforward to perform and will be 
reported elsewhere. 

2.5. Interpolation on the sphere 

In practice we have used the HEALPix 9 (Gorski et al. 
2005) scheme to represent fields on the sphere. At various 
stages of the lensing calculation, an accurate algorithm 
for interpolation on the sphere becomes a necessity. In 
the effective lensing approximation, the original positions 
of the rays will in general be off pixel centers. This im- 
plies that the lensed CMB field is essentially generated by 
sampling the unlensed CMB surface at points which are 
usually not pixel centers. Hence, obtaining the lensed 
CMB field is essentially an interpolation operation. In 
the case of multiple lensing planes, it is again obvious 
that (except for the first plane, on which we can shoot 
rays at pixel centers by way of convenience) the deflection 
field itself has to be evaluated at off-center points on all 
subsequent planes. So, together with the interpolation 
of the temperature map. we need to go between spin-0 
and spin-1 fields on an arbitrary grid. Therefore, one 
needs, in general, a spherical harmonic transform algo- 
rithm that can deal with an irregular grid on the sphere. 

For this purpose, we adopt the Non-isolatitude Spheri- 
cal Harmonic Transform (NISHT) algorithm proposed by 
Hirata et al. (2004); details of the algorithm can be found 
in Appendix A of that paper. Here we have reproduced 
the key equations for clarity, and described the salient 
features of the general algorithm with special attention 
to aspects which are relevant for the current application. 

The basic operation for generating the lensed CMB 
maps can be broken up into two steps: 

LI. generating the deflection field on the sphere at 
points where the rays land from the previous plane, and 

L2. sampling the unlensed CMB surface at the source- 
plane positions of the rays to generate the lensed CMB 
field. 

Of course, in case of the effective lensing simulation, 
one can conveniently generate the deflection field at the 
pixel centers in step LI above. As step L2 is a series of 
operations involving scalars and therefore conceptually 
simpler, we shall explain the NISHT algorithm in relation 
to this step. Step LI, which involves spin-1 fields on the 
sphere, is conceptually similar to the spin-0 case. 

The problem in step L2 is that we know the CMB tem- 
perature field T(n) on the HEALPix grid {n} as well as 
the source-plane positions of the rays {n'} on the po- 
lar cap, and we want to sample the CMB field at {n'}. 
Suppose, by applying the steps for spherical harmonic 
analysis (as will be described later), we have the spheri- 
cal harmonic components, Tg m of the temperature field. 
Now, we need to synthesize the field using these T^ m 's at 
the points {n'}. This operation can be formally written 
as 



= E E T tm Y em (n'), 



(29) 



(=0 m= 



where £ max 

resolution of the HEALPix grid as 



is the Nyquist multipole and is set by the 

ax — ?r/ \Z^pix 



9 http://healpix.jpl.nasa.gov 
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(cf. equation 40). This synthesis operation can be split 
into the following four steps (Eqns. 30 through 35 are 
essentially reproduced for completeness from Appendix 
A of Hirata et al. 2004): 

1. Coarse Grid Latitude Transform 

As the hrst step, we perform a transform in the lat- 
itude direction on an equally spaced set of points, 
(6 = ■na/L,<j) — 0), where a is an integer in the 
range < a < L and L is a small integral multiple 
of some power of two such that L > £ max : 



T m (8 = — 71") = E TtmYir, 

e=\ m \ 



-7T,0 = O). (30) 



The above calculation involves 0(€ max L) opera- 
tions. 

2. Refinement of Latitude Grid 

In this step we reduce the 9 grid spacing from ^7r 
to jj-k where L' > L. We take advantage of the 
fact the sampling theorem can be applied to a lin- 
ear combination of spherical harmonics which is 
band limited (£ < ^ max ) in the multipole space, 
and hence can be written as a Fourier sum, 



,(*)= E C m , n e 



in9 



(31) 



We determine the coefficients C TOi „ via a fast 
Fourier Transform (FFT) of length 2L and eval- 
uate T m {9 = jj-k) using an inverse FFT of length 
2L 1 . This step saves us the expensive generation of 
Associated Legendre Polynomials on the finer grid. 
Each FFT requires 0(£ max L log(L)) operations. 

3. Projection onto Equicylindrical Grid 

Next, we perform the standard SHT step of taking 
an FFT in the longitudinal direction to generate 
T(0 fa), 



= 77^) 



E T ^ e 



(32) 



After this step we have synthesized the map on an 
Equicylindrical projection (ECP) grid. The opera- 
tion count for this step is 0(L' 2 log L') and the total 
operation count including this step is 0(^ max 3 ). 

4. Interpolation onto the final grid 

In the final step, given a required position n', we 
find the nearest grid point in the ECP grid and de- 
termine the fractional offset, {8a, Sj) between the 
two points, 

a + Sa = L' — ; 7 + 67 = L' — . (33) 

7T 7T 

Then we perform a two dimensional polynomial in- 
terpolation using (2K) 2 points around the nearest 
grid point, obtaining the value at the required point 



K 



K 
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a + [i 7 + v 



v=-K+l 



(34) 



with the weights computed using Lagrange's inter- 
polation formula, 



w p (5) 



(-l) K - p 



K 



(K-p)\(K-l + P )\{5-p) 



n 



-A'+l 



(35) 



The inverse of the synthesis operation described above 
is the analysis operation, in which the spherical har- 
monic coefficients of a map defined on an irregular grid 
is needed. This can be thought of as the transpose of 
the above operations applied in reverse, and hence can 
be accomplished in an equal number of steps. 

The above algorithm can be easily extended to deal 
with vector and tensor fields on the sphere. For a vec- 
tor (spin 1) field, the natural basis of expansion are the 
vector spherical harmonics, 



■yV „. 

x e m — 



1 



■ Ira " 



V^TT) 
1 



n x VYfr, 



(36) 



where the superscripts V and A represent the "vector- 
like" and the "axial-vector-like" components, respec- 
tively. In terms of these a vector field v(n) can be ex- 
panded as, 

v ( fi ) = EE ^» Y L(n) + AtmYLW- (37) 

e=o m=-l 

Therefore, given the (Ve m ,Ae m ) components one can go 
through the analogs of the above steps for the scalar field 
synthesis. In fact, to accomplish step LI of the lensing 
algorithm, we go from the convergence field K to the 
deflection field, 



a. 



e=o m=~ 

£max i 

=-E E 



(■in ■in 



e=o m=-i 



£(£ + 1) 



=-E E 



£=0 m-- 



■yV 
i 1 lm- 



(38) 



Therefore, we go from the K field on the polar cap to the 
spherical harmonic components Kp m using the analysis 
algorithm for scalar fields; then we divide the result by 
\J l{t + l)/2. This defines the vector field harmonic com- 
ponents as {Ve m ,A( m ) = (-2K im / y / £(£+ 1),0) from 
which we synthesize the deflection field at the required 
points. 

The accuracy of the interpolation can be controlled by 
two parameters: the rate at which the finer grid over- 
samples the field i.e. the ratio L'/^ max , and the order 
of the polynomial K used for the interpolation. Increas- 
ing either of these increases the accuracy. In this paper 
we have used L' = 4^ max and K = 10, which yields 
a fractional interpolation accuracy per Fourier mode of 

~io- 9 . 
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3. GENERATION OF THE LENSING PLANES 

An Af-body dark matter simulation was performed to 
generate the large scale structure; this same simulation 
has been discussed in Schgal et al. (2007) and Bode 
et al. (2007), so we refer the reader to these papers 
for more details. Briefly, a spatially flat ACDM cos- 
mology was used, with a total matter density parame- 
ter fi m = 0.26 and vacuum energy density ft a = 0.74. 
The scalar spectral index of the primordial power spec- 
trum was set to n s — 0.95 and the linear amplitude 
normalized to cr 8 = 0.77. The present day value of 
the Hubble parameter H = 72 kms _1 Mpc _1 . A pe- 
riodic box of size L = 1000 /i _1 Mpc was used with 
N = 1024 3 particles; therefore the particle mass was 
m p = 6.72 x 10 10 /i _1 M Q . The cubic spline softening 
length was set to 16.28 /i _1 Mpc. 

3.1. From the box to the sphere 

We create the lensing planes on-the-fly from the TV- 
body simulation. At each large time step (set by a 
Courant condition such that no particle moves more than 
~ 122/i~ 1 kpc in this time) the positions and velocities of 
the particles in a thin shell are saved. The mean radius 
of the shell is the comoving distance to the redshift at 
that time, and the width (a few ft _1 Mpc) corresponds to 
the time step. Each shell is centered on the origin of the 
simulation and covers one octant of the sky (x, y, z > 0). 
Note that for shells with radii greater than the simula- 
tion box size, periodic copies of the box are used. Thus 
a given structure will appear more than once in the full 
light cone, albeit at different times and viewed from dif- 
ferent angles. We then Eulcr rotate the coordinate axes 
so that the new z-axis passes through the centroid of the 
octant. This is done to make the centroid correspond 
to the North Pole on the HEALPix sphere. We use the 
HEALPix routine vec2pix to find the pixel that con- 
tains the particle's position on the sky. We then place 
the mass of the particle into that pixel by assigning to 
it the surface mass density E p = m p /Q p i x , where f2 p i x is 
the area of a pixel in steradians (cf. equation 40). Thus, 
if n particles fall inside the beam defined by a pixel, then 
the pixel ends up having a surface mass density of nT, p . 
To simplify the geometry, we save only those particles 
which fall inside a Polar Cap like region defined by the 
disc of maximum radius that can be cut out from the 
octant (see Fig. 3). 

By the end of the run, 449 such planes were produced 
from the simulation, spanning z = 4.0 to z = 0. As these 
are far too many planes for the purpose of lensing, we 
reduce them into ~ 50 planes by dividing up the original 
planes into roughly equal comoving distance bins and 
adding up the surface mass density pixel by pixel for all 
planes that fall inside a bin to yield a single plane per 
bin. Hereafter, we shall refer to the original planes from 
the TPM run as the TPM-planes and the small number 
of planes constructed by projecting them as the lensing- 
planes. 

The angular radius of the Polar Cap is given by # cap = 
arccos(2/v / 6), and the solid angle subtended by it is 



2tt(1 - cos6» ca p) = 1.981 sr = 3785 sq - dcg. Due 



to pixelation, the true total area f2 cap of the iV ca p pixels 

9 http://www.jportsmouth.com/code/CMBview/cmbview.html 




Fig. 3. — Illustration of the Polar Cap geometry. The 
figure shows a 3-D rendering of the sphere using the CMB- 
VIEW a software, looking down towards the North Pole. The lightly 
shaded triangular region correspond to the positions of the parti- 
cles in the octant from the iV-body simulation box at a typical 
time step. The darker dots define what we call the Polar Cap in 
the text. The surface mass density of the pixels inside this Polar 
Cap region are saved in shells out to 2=4. 

that make up the Polar Cap is not exactly equal to f2 C ap, 
but rather 

(39) 

We will denote the surface mass density in pixel p as a p 
which has units of mass per steradian. 

In HEALPix , the resolution is controlled by the pa- 
rameter NSIDE, which determines the number, -/V p ; x of 
equal area pixels into which the entire sphere is pixe- 
lated, through the relation N pix = 12 x NSIDE 2 , so that 
the area of each pixel becomes, 



^cap — -^cap ^pix • 



n 



pix 



4-7T 
1 » nix 



steradians. 



(40) 



The angular resolution is often expressed through the 
number 8 ICS = y/Qpix . It is also useful to define the 
fraction area of the sphere covered by the polar cap as, 



sky 



cap 



47T 



(41) 



For results presented in this paper the resolution pa- 
rameter NSIDE was set to 4096 , which corresponds to 
an angular resolution of 0/896 . 

3.2. From surface density to convergence 

To construct the quantities required for lensing, we 
first convert the surface mass density maps into surface 
over-density maps AE as defined in equation (14). It is 
straightforward to obtain the if-maps defined in equa- 
tion (15) from the above map. Finally, equation (20) is 
used to obtain the effective convergence map on the Polar 
Cap. It is evident that the convergence map constructed 
out of the simulated lensing planes in this way will only 
contain the contribution from large scale structure up 
to the redshift of the farthest lensing plane (z = 4.05). 
However, to accurately lens the CMB we need to add 
in the contribution from higher redshifts up to the last 
scattering surface. We do this by generating a Gaussian 
random field from a theoretical power spectrum of the 
convergence between z = 4.05 and z = z CMB , computed 
from the matter power spectrum obtained using CAMB, 
and adding it onto the convergence map from the TPM 
simulation. 
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3.3. The unlensed CMB map 

We used the synfast facility in HEALPix to gen- 
erate the unlensed CMB map. This takes as an in- 
put a theoretical unlensed power spectrum and synthe- 
sizes a Gaussian-random realization of the unlensed CMB 
field. For computing the theoretical power spectrum we 
have used the publicly available Boltzmann transfer code 
CAMB 10 , with the same set of cosmological parameters 
as used for the large scale structure simulation. 

4. MEASURING ANGULAR POWER SPECTRA 

At several stages we compute the power spectra of the 
maps to compare with theory. For example, to verify 
that we have created the convergence map correctly, the 
angular power spectrum of the k map is computed and 
compared to the theory. Also, we do the same for the 
lensed map on the polar cap. We use the map2alm facil- 
ity of HEALPix to perform a spherical harmonic decom- 
position of a map T(h) on the Polar Cap. The result- 
ing spherical harmonic components, i.e. T^'s, are then 
combined to obtain the pseudo-power spectrum, 



1 



21+1 ^ 



\ 1 im 



(42) 



There are two effects that need to be taken into ac- 
count before comparing the above result with theory, 
namely the finite pixel size, signified by the superscript, 
"pix" and the incomplete sky coverage, represented by 
the tilde. 

To simplify the following discussion of pixelation ef- 
fects, for the moment we shall ignore the effect of the in- 
complete sky coverage. Also, we shall use the shorthand 
notation T,g m to denote the sum J2eLo^2m=~e- Due to 
the finite pixel size, a field realized on the HEALPix 
sphere is a smoothed version of the true underlying field, 
i.e. the value of the field in pixel i is given by 



j>pix(i) 



d 2 nw^(h)T(h), 



(43) 



where u;W is the window function of the i-th pixel as is 
given by 



fi pix inside pixel i, 
elsewhere. 



(44) 



Expanding the true field T in terms of spherical harmon- 
ics as 



T(n) =J2 T emY em (n), 



im 



we have 



where 



tm 



d 2 hw^(h)Y lm (h) 



(45) 



(46) 



is the spherical harmonic transform of the pixel window 
function. In the HEALPix scheme, due to the azimuthal 
variation of the pixel shapes over the sky, especially in 

10 http://www.camb.info 



the polar cap area, a complete analysis would require 
the computation of these coefficients for each and every 
pixel. However, even for a moderate NSIDE, this calcula- 
tion becomes computationally unfeasible. Therefore, it is 
customary to ignore the azimuthal variation and rewrite 
equation (46) as 



(47) 



where one defines an azimuthally averaged window func- 
tion 

P 1 V2 



(*) 
W) = 



47T 



(2£+l) 



E 



(48) 



From equations (47) and (45) it immediately follows that 
the estimate of the power spectrum of the pixelated field 
is given by 

Cf x = w 2 (T em T e * m ) = w 2 C e (49) 
where one defines the pixel averaged window function, 

1/2 

{wff ) . (50) 



f 



we 



AW 



E 

i=0 



This function is available for I < 4 x NSIDE in the 
HEALPix distribution. We take out the effect of the 
pixel window by dividing the computed power spectrum 
by the square of the above function. Coming back to 
the case at hand, where we have both pixel and incom- 
plete sky effects, we recover the power spectrum Ce after 
correcting for the pixel window function in this manner. 

The second and more important effect that one needs 
to take into account results from the fact that our field 
is defined only inside the polar cap. This is equivalent to 
multiplying a full sky map with a mask which has value 
unity inside the polar cap and zero outside. As is well 
known, such a mask leads to a coupling between various 
multipoles, leading to a power spectrum which is biased 
away from the true value. As this effect tends to move 
power across multipoles, the problem is more acute for 
highly colored power spectra like the CMB. 

Let us denote the effective all-sky mask with W, where 



W(n) = {l 



inside the polar cap, 



elsewhere. 



(51) 



The spherical harmonic components of the masked field 
is therefore given by 

t lm = J d 2 hT(n)W(n)Y; m (h) (52) 
■■Y,Tt'm' I d 2 hY llm ,(n)W(h)Y; m (n) (53) 

I'm' J 



and the measured power spectrum by (see for example 
Hivon et al. 2002) 



^ - (24 + l) E (^"H^r 
mi=-<i 



(54) 
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Fig. 4. — Effect of apodization of the window function. The con- 
tinuous line is the unlensed CMB power spectrum and the dashed 
line is the lensed one. Both have been scaled by , the effective 
fractional sky coverage (see text). The gray filled and open circles 
labeled "Tophat", represent, respectively, the theoretical unlensed 
and lensed power spectra convolved with a window function that 
is unity inside the polar cap and zero outside. The black filled 
and open circles represent the same quantities, but in the case of 
a window which is apodized at the edge of the polar cap, as dis- 
cussed in the text. Aliasing of power to higher multipoles due to 
mode coupling is significantly reduced in the latter case. We use 
the apodized window to mask the polar cap maps for computing 
various power spectra, and use the corresponding theory power 
spectrum convolved with the same window for comparing our re- 
sults with theory. {Inset: The mode coupling matrix as a 
function of l\ , for £2 = 3000, showing the reduction in the power 
in off-diagonal elements as a result of apodization.) 

where Ce 2 is the true power spectrum and M is the mode 
coupling matrix given by 

.^Mq^pJJ oV (55) 

with the power spectrum of the mask defined as 

^ = ^TI^ |v ^ m|2 ' (56) 

m 

Wim being the spherical harmonic components of the 
mask W(n). 

For a polar cap with angular radius 9, this function is 
analytically known to be (Dahlcn & Simons 2007) 

Wf P - {2e l 1)2 [P/-i(cos0) - P, +1 (cos9)] 2 . (57) 

where Pg is a Lcgendre Polynomial of order £ and 

The window function in equation (51) corresponding to 
the polar cap is a "tophat" in the sense that it abruptly 
falls to zero at the edge. The power spectrum (equation 
(57)) of this window has an oscillatory behavior showing 
a lot of power over a large range of multipoles, an effect 
sometimes called ringing. Ringing causes the mode cou- 
pling matrix, Mu> to develop large off-diagonal terms, 
as illustrated in Fig. 4, and consequently the value of 



the measured power spectrum at any multipolc (equation 
(54)) has non-trivial contributions from many neighbor- 
ing multipoles. This causes the measured power spec- 
trum to be biased, and its effect is particularly evident 
for power spectra with a sharp fall-off such as the CMB. 
As is evident from Fig. 4, the effect of mode coupling due 
to the polar cap becomes a serious problem for the lensed 
and unlensed power spectra starting at moderately low 
multipoles (£ <~ 2000). Although in principle one could 
compare the measured power spectrum with a theoretical 
power spectrum which has been convolved with the same 
mode coupling matrix, the effect is so strong in this case 
that the lensed and unlensed spectra almost overlap each 
other. This problem can be mitigated in principle by in- 
verting a binned version of the mode coupling matrix 
and thereby decorrelating the power spectra. However 
an easier and less computationally expensive solution can 
be achieved in the following manner. 

The off diagonal terms of the mode coupling matrix 
can be reduced significantly by apodizing the polar cap 
window function. Parenthetically, we note that there ex- 
ists a general method of generating tapers on a cut-sky 
map, so as to minimize the effect of mode coupling. This 
is referred to as the multi-taper method (S. Das, A. Ha- 
jian and D. N. Spergel, 2007, in preparation; Dahlen & 
Simons 2007). However, for our purpose, it suffices to 
define a simpler apodizing window as 

fi for e<e <e C!ip 

W(n) = \ sin(f §^§-J for 9 < 9 < 9 cap (58) 
I for 9 > 9 cap . 

The power spectrum of this window can be easily com- 
puted using HEALPix , and thus the mode coupling ma- 
trix can be readily generated using equation (55). We 
found that an apodization window with (# cap — 9q) — 1.2 
degree, corresponding to <~ 80 pixels, works extremely 
well without eating into too much of the map. A sec- 
tion of the mode coupling matrix and the correspond- 
ing convolved power spectrum are displayed in Fig. 4. 
This figure shows that the power spectrum convolved 
with the apodized window function has negligible mode 
coupling. Parenthetically, it is interesting to note that 
simply scaling the theory power spectrum by the frac- 
tion of the sky covered, / s k y , seems to do a good job in 
mimicking the effect of the partial sky coverage, at least 
for the lower multipoles. In fact, this approximation is 
an exact result for a white power spectrum. However, 
when the window is apodized, the effective area of cov- 
erage, = J W 2 (h)d 2 h/4:'K, goes down a little (by 

<~ 2.5% for our apodization). We use f^ y scaled theory 
power spectra only in some plots in this paper. For the 
analysis, we perform the full mode coupling calculation. 
Therefore, when comparing the power spectrum of some 
quantity defined on the polar cap with theoretical predic- 
tions, we first multiply the map by the apodized window 
and compare the resulting power spectrum with the the- 
oretical power spectrum mode-coupled through the same 
weighting window. 

5. RESULTS 

We illustrate the algorithm with an effective lensing 
simulation at the HEALPix resolution of NSIDE = 
4096. Since some rays end up outside the polar cap after 
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Fig. 5. — The Polar Cap map obtained after subtracting the 
unlensed CMB map from the lensed CMB map. To enhance the 
contrast, we have remapped the color scale to the range (— 2cr, 2a), 
a being the standard deviation of the map. 

lensing, we have actually used an unlensed CMB real- 
ization (using the synf ast facility of HEALPix ) on an 
area larger than our fiducial polar cap to accommodate 
those rays. As the gradient of the lensing potential is 
ill defined at the edge of the polar cap, we ignore a ring 
of pixels near the edge of the lensed map for all subse- 
quent analyses. It is particularly instructive to look at 
the difference of the lensed and the unlensed maps, as 
shown in Fig. 5, as it shows the large scale correlations 
imprinted on the CMB due to the large scale modes in 
the deflection field. 

We compute the angular power spectrum of the lensed 
and unlensed CMB maps using an apodized weighting 
scheme as discussed in §4. The resulting power spec- 
tra are displayed in Fig. 6 for the entire range of multi- 
poles analyzed, and are compared with the mode-coupled 
theoretical power spectra. The theoretical lensed CMB 
power spectrum used for the calculation was generated 
with the CAMB code, using the all-sky correlation func- 
tion technique (Challinor & Lewis 2005) and including 
nonlinear corrections to the matter power spectrum. In 
Fig. 7, we show a zoomed-in version of the lensed power 
spectrum, in the multipole range 500 < I < 3500. From 
visual inspection of these plots it is evident that the sim- 
ulation does a good job in reproducing the theoretically 
expected lensed power spectrum, at least in the range 
of multipoles over which the computation of the theo- 
retical power spectrum is robust. We defer a detailed 
comparison of the simulation to the theory to §6.3. 

6. TESTS 

6.1. Tests for the mass sheets 

In this section we perform some sanity checks to ensure 
that the projection from the simulation box onto the Po- 
lar Cap has been properly performed. We first test that 



Fig. 6. — The lensed and unlensed CMB angular power spectra 
obtained from the simulation compared with the theoretical mod- 
els. The red and orange dots represent, respectively, the lensed and 
unlensed angular power spectra obtained from the polar cap using 
the methods described in §4. The solid black curve signifies the 
theoretical unlensed power spectrum taking into account the mode 
coupling due to the apodized polar cap window function. The blue 
solid curve represents the same for the lensed power spectrum. 
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Fig. 7. — Lensed CMB angular power spectrum in the multipole 
range 500 < £ < 3500 obtained from the simulation compared with 
the theoretical model. The red dots represent the lensed angular 
power spectrum obtained from the polar cap using the methods de- 
scribed in §4. The solid black curve signifies the theoretical lensed 
CMB power spectrum taking into account the mode coupling due 
to the apodized polar cap window function. The dotted black curve 
represents the same for the theoretical unlensed power spectrum 
and is shown here for contrast to the lensed case. 



the total mass in each slice is equal to the theoretical 
mass expected from the mean cosmology, the later being 
given by 

M^cr = fWcritfttcapA?? (59) 
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Fig. 8. — The mass in the lensing-sliccs compared with that 
expected from theory. 

where A77 is the comoving thickness of the slice at a co- 
moving distance fj. We compare this quantity with 

Neap 

M sVlcc = ^ o-jfipix (60) 

i=i 

which is the total mass on the plane from the simulation. 
The percentage difference between the two is depicted in 
Fig. 8 for the lensing-planes. Notice that the agreement is 
good to within 0.5% for the high redshift planes, in which 
the solid angle O cap corresponds to a large comoving area. 
For low redshifts there are large variations due to the fact 
that matter is highly clustered and £! cap corresponds to a 
small comoving area. These fluctuations at low redshift 
represent the chance inclusion or exclusion of large dark 
matter halos within the light cone. 

Next, we make sure that the probability density func- 
tion (PDF) of the surface mass density is well behaved 
for each plane, and is well modeled by analytic PDFs 
such as the lognormal (Kayo et al. 2001; Taruya et al. 
2002) or the model proposed by Das & Ostriker (2006). 
In Fig. 9 we show these two models over-plotted on the 
PDFs drawn from the forty-five lensing-planes. 

The model of Das & Ostriker (2006) is a better fit 
to the simulation than the lognormal, especially at high 
surface mass density. Note in that paper the authors used 
the first year WMAP parameters, whereas the present 
simulation is run with the WMAP 3-year parameters, 
including a significantly different er 8 . The fact that the 
model still represents the simulation well suggests that 
it is quite general. 

6.2. Tests for the convergence plane 

As described in §3.2, the effective convergence plane 
was produced by a two step process. First, we computed 
the effective convergence plane by weighting the surface 
mass density planes from the simulation with appropri- 
ate geometrical factors. Let us call it the map Mi. This 
map, therefore, includes contribution from the large scale 
structure only out to the redshift of the farthest TPM 
plane, z — 4.05. Next we added in the contribution from 
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Fig. 9. — The probability density function (PDF) of the surface 
mass density in the lensing-planes (circles) compared with the log- 
normal (dashed line) and the Das & Ostriker (2006) model (solid 
line). 

z > 4.05, by generating a Gaussian random realization of 
the effective convergence from a theoretical power spec- 
trum (the map M2). Therefore the final convergence 
map is simply (Mi + M 2 ). It is interesting to compare 
the power spectrum of the map Mi with that expected 
from theoretical considerations. Since CMB lensing is 
most sensitive to large scale modes, we should make sure 
that these modes were realized correctly in our simulated 
convergence plane. Incidentally, these scales are also lin- 
ear to mildly nonlinear. Therefore, we should expect 
the power spectrum of the convergence map to be well 
replicated by the theoretical prediction at least in the 
quasi-linear range of multipoles (t < 2000) where simple 
non-linear prescriptions suffice. 

In order to compute the theoretical power spectra for 
the maps Mi and (Mi + M 2 ), we used the Limber 
approximation to project the matter power spectrum 
P(fc, 77) computed from CAMB. The Limber approxi- 
mation simplifies the full curved-sky calculation, and is 
valid for I > 10. Since for lower values of the multi- 
pole we have few realizations of the convergence modes, 
the power spectrum computed from the simulated map 
is noisy in this regime, rendering it practically useless for 
comparison with theory. Therefore, an accurate compu- 
tation of the theoretical convergence power spectrum for 
these lowest multipoles is unnecessary, and the Limber 
calculation suffices. Under the same approximation, the 
shot noise contribution to the convergence field can be 
computed as 

-shot \^ A „ ( 3 o ti , ^ ('ns-Vj)Vj H$\ 2 1 

(61) 

where rij — Nj/(r]jAr]jSl ca p), Nj being the total number 
of particles in the j-th shell. 

We compute both a linear and a non-linear version 
of the convergence power spectrum, where the latter in- 
cludes non-linear corrections to the matter power spec- 
trum from a halo model based fitting formula (Smith 
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Fig. 10. — Power spectrum of the effective convergence map Mi 
produced from the simulated lensing planes alone. The red line 
shows the power spectrum computed from the convergence map 
and the black solid line represents the theoretical power spectrum 
with non-linear corrections. The power spectrum is corrected for 
the shot noise contribution (see text) which is displayed as the 
dotted line. The black dashed line corresponds to the linear theory 
power spectrum. All theory power spectra are mode-coupled with 
the apodizing window. 
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Fig. 11. — Power spectrum of the effective convergence map 
(Mi +M2) after adding in high redshift contribution. The red line 
shows the power spectrum computed from the convergence map 
and the black solid line represents the theoretical power spectrum 
with non-linear corrections. The power spectrum is corrected for 
the shot noise contribution (see text) which is displayed as the 
dotted line. The black dashed line corresponds to the linear theory 
power spectrum. All theory power spectra are mode-coupled with 
the apodizing window. 

et al. 2003). We plot the power spectrum computed 
from the simulated convergence map Mi, and the cor- 
responding theoretical power spectra, in Fig. 10. As is 
evident from the figure, the simulated power spectrum 
is in accord with the linear theory power spectrum for 
£ < 300, beyond which the effect non-linearities creep in. 
However, it is impressive that the non-linear corrections 



to the power spectrum are in good agreement with the 
simulation up to relatively high multipoles. The same 
quantities are plotted for the convergence map out to 
the redshift of the CMB in Fig. 11. We find in both 
cases that beyond multipoles of ~ 6000 the simulation 
contains more non-linear power than predicted by the 
theory. 

6.3. Tests for the lensed CMB map 

Since CMB lensing is essentially a remapping of points, 
the one-point statistics should remain unaffected by the 
lensing. We check for this by drawing up the one-point 
PDF's of the unlensed and lensed maps, and find them 
to be consistent to within 0.8%. Next, we compare the 
power spectrum of the simulated lensed map (cf. Fig- 
ures 6 & 7) with the theoretical predictions as computed 
with the CAMB code. For a quantitative comparison, 
we consider a range of multipoles (500 < I < 3500) in 
the acoustic regime. We do not consider the lower multi- 
poles as they exhibit negligible lensing effect. We found 
that for a fixed input cosmology, the tail (£ > 3500) of 
the lensed CMB power spectrum predicted by CAMB de- 
pends somewhat sensitively on input parameters, specif- 
ically the combination krj mstX: which controls the max- 
imum value of the wavenumber for which the matter 
power spectrum in computed. However, the lensed power 
spectrum from CAMB is robust towards changes in the 
auxiliary input parameters for the range of multipoles, 
I < 3500. Also, the lensed CMB multipoles beyond this 
range couple to relatively small scale modes of the deflec- 
tion field where our simulation has more power than ex- 
pected from non- linear theory. In fact, beyond I ~ 4000 
the simulated power spectrum is found to deviate sys- 
tematically from the theoretical spectrum. 

As the simulated power spectrum, C| lm , was computed 
using an apodized window as described in §4, the appro- 
priate theoretical curve to compare this result with is the 
power spectrum from CAMB after it has been convolved 
with the coupling matrix defined by the same weighting 
scheme (cf. equation 54 ) . We denote the latter quantity 

by C[ hcoTy . In order to facilitate the comparison we bin 
the raw spectrum in I. In the multipole range consid- 
ered (500 < £ < 4000), the quantity C e = £ 4 C e is flat 
(see Fig. 7) and therefore a better candidate for binning. 
We denote the difference between the simulated and the 
theoretical version of this quantity by 



SCi = c? im - c, 



theory 
I 



(62) 



For each of the N bins indexed by 6, we compute the 
mean, SCb, and the sample variance, s^, of the obser- 
vations falling inside that bin. In order to account for 
that fact that cosmic variance errors will be higher in 
our case due to incomplete sky coverage, we define an 
effective variance as a\ — s\/ 'f^y- 

We quantify the goodness of fit between the simulation 
and the model by defining a \ 2 statistic as 



x 2 = 



N 

E 

b=l 



(5C b f 



(63) 



We perform the \ 2 analysis by uniformly binning the 
power spectra in the range 500 < £ < 4000 into 52 bins 
with a bin width of A£ — 60. The binned values along 
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Fig. 12. — Difference between the simulated and the theoretical 
binned power spectrum for lensed CMB. 

with the error bars are displayed in Fig. 12. We find a 
value x 2 = 52.93, suggesting an appreciable agreement 
with the theory. 

7. CONCLUSIONS 

In this paper, we have put forward an algorithm for 
end-to-end simulation of the gravitational lensing of the 
cosmic microwave background, starting with an TV-body 
simulation and fully taking into account the curvature of 
the sky. The method is applicable to maps of any ge- 
ometry on the surface of the sphere, including the whole 
sky. Our algorithm includes prescriptions for generat- 
ing spherical convergence planes from an N-body light 
cone and subsequently ray-tracing through the planes to 
simulate lensing. The central feature of the algorithm is 
the use of a highly accurate interpolation method that 
enables sampling of both the deflection fields on inter- 
mediate lensing planes and the unlensed CMB map on 
an irregular grid. We have provided a detailed descrip- 
tion of both a multiple plane ray tracing and an effec- 
tive lensing version of the algorithm. The latter setting 
has been used to illustrate the algorithm, by generat- 
ing an ~1' resolution lensed CMB map. We have com- 
pared the power spectra of the effective convergence map 
and the lensed CMB map with theoretical predictions, 
and have obtained good agreement. After this paper 
was completed, Fosalba et al. (2007) described a similar 
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method of producing convergence maps in spherical ge- 
ometry, and Carbone et al. (2007) also described their 
techniques for simulating CMB lensing maps. The latter 
used broadly similar techniques to those described here, 
although they used a different method to obtain the de- 
flection field. 

Applications of the algorithm can be manifold. The 
associated large scale structure planes can be populated 
with tracers of mass and foreground sources, in order to 
simulate cross-correlation studies and to investigate the 
effects of contamination. This lensing portion of the al- 
gorithm can be applied to generate lensed maps in large 
scale structure simulations that produce spherical shells 
(Fosalba et al. 2007). The multiple plane algorithm can 
be particularly useful, after trivial modifications, in sim- 
ulating weak lensing of galaxies or the 21-cm background 
on a large sky. The lensed CMB maps can be used as 
inputs to telescope simulators for projects such as ACT 
and PLANCK, and will help in the analysis and inter- 
pretation of data. We intend to release all-sky high res- 
olution lensed CMB maps made using this algorithm in 
near future. 



SD would sincerely like to thank his advisor, David 
Spergel for suggesting the key ideas of the project and 
for his continuous guidance and encouragement through- 
out its development. SD is specially grateful to Chris 
Hirata for generously providing the NISHT code and for 
numerous useful discussions. We thank Joanna Dunk- 
ley for careful reading of the manuscript and thoughtful 
suggestions. We would also like to thank the referee for 
his thoughtful suggestions. SD acknowledges the sup- 
port from the NASA Theory Program NNG04GK55G 
and the NSF grant AST-0707731. This research was fa- 
cilitated by allocations of advanced computing resources 
from the Pittsburgh Supcrcomputing Center and the Na- 
tional Center for Supercomputing Applications. In addi- 
tion, computational facilities at Princeton supported by 
NSF grant AST-0216105 were used, as well as high per- 
formance computational facilities supported by Prince- 
ton University under the auspices of the Princeton Insti- 
tute for Computational Science and Engineering (PIC- 
SciE) and the Office of Information Technology (OIT). 
Some of the results in this paper have been derived using 
the HEALPix (Gorski et al. 2005) package. 



Hirata, C. M., Padmanabhan, N., Seljak, U., Schlegel, D., & 

Brinkmann, J. 2004, Phys. Rev. D, 70, 103501 
Hirata, C. M. & Seljak, U. 2003, Phys. Rev. D, 67, 043001 
Hivon, E., Gorski, K. M., Netterfield, C. B., Crill, B. P., Prunet, 

S., & Hansen, F. 2002, ApJ, 567, 2 
Hu, W. 2000, Phys. Rev. D, 62, 043007 
Hu, W., Holz, D. E., & Vale, C. 2007, ArXiv e-prints, 708 
Hu, W., Huterer, D., & Smith, K. M. 2006, ApJ, 650, L13 
Hu, W. & Okamoto, T. 2002, ApJ, 574, 566 
Jain, B., Seljak, U., & White, S. 2000, ApJ, 530, 547 
Kayo, I., Taruya, A., & Suto, Y. 2001, ApJ, 561, 22 
Komatsu, E., Kogut, A., Nolta, M. R., Bennett, C. L., Halpern, 

M., Hinshaw, G., Jarosik, N., Limon, M., Meyer, S. S., Page, L., 

Spergel, D. N., Tucker, G. S., Verde, L., Wollack, E., & Wright, 

E. L. 2003, ApJS, 148, 119 
Lcsgourgues, J., Perotto, L., Pastor, S., & Piat, M. 2006, 

Phys. Rev. D, 73, 045021 
Lewis, A. 2005, Phys. Rev. D, 71, 083008 



14 



Das & Bode 



Lewis, A. & Challinor, A. 2006, Phys. Rep., 429, 1 
Marian, L. & Bernstein, G. M. 2007, ArXiv e-prints, 710 
Okamoto, T. & Hu, W. 2003, Phys. Rev. D, 67, 083002 
Peiris, H. V. & Spergel, D. N. 2000, ApJ, 540, 605 
Sehgal, N., Trac, H., Huffenberger, K., & Bode, P. 2007, ApJ, 664, 
149 

Seljak, U. 1996, ApJ, 463, 1 

Smith, K. M., Hu, W., & Kaplinghat, M. 2006a, Phys. Rev. D, 74, 
123002 

Smith, K. M., Zahn, O., & Dore, O. 2007, Phys. Rev. D, 76, 043510 
Smith, R. E., Peacock, J. A., Jenkins, A., White, S. D. M., Frenk, 

C. S., Pearce, F. R., Thomas, P. A., Efstathiou, G., & Couchman, 

H. M. P. 2003, MNRAS, 341, 1311 
Smith, S., Challinor, A., & Rocha, G. 2006b, Phys. Rev. D, 73, 

023517 

Spergel, D. N., Bean, R., Dore, O., Nolta, M. R., Bennett, C. L., 
Dunkley, J., Hinshaw, G., Jarosik, N., Komatsu, E., Page, L., 
Peiris, H. V., Verde, L., Halpern, M., Hill, R. S., Kogut, A., 
Limon, M., Meyer, S. S., Odegard, N., Tucker, G. S., Weiland, 
J. L., Wollack, E., & Wright, E. L. 2007, ApJS, 170, 377 



Spergel, D. N., Verde, L., Peiris, H. V., Komatsu, E., Nolta, M. R., 
Bennett, C. L., Halpern, M., Hinshaw, G., Jarosik, N., Kogut, 
A., Limon, M., Meyer, S. S., Page, L., Tucker, G. S., Weiland, 
J. L., Wollack, E., & Wright, E. L. 2003, ApJS, 148, 175 

Taruya, A., Takada, M., Hamana, T., Kayo, I., & Futamase, T. 
2002, ApJ, 571, 638 

Vale, C, & White, M. 2003, ApJ, 592, 699 

Zaldarriaga, M. 2000, Phys. Rev. D, 62, 063510 

Zaldarriaga, M. & Seljak, U. 1999, Phys. Rev. D, 59, 123507 



